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ABSTRACT 

We introduce a lattice Boltzmann computational scheme capable of modeling ther- 
mohydrodynamic flows of monatomic gases. The parallel nature of this approach 
provides a numerically efficient alternative to traditional methods of computational 
fluid dynamics. The scheme uses a small number of discrete velocity states and a 
linear, single-time-relaxation collision operator. Numerical simulations in two dimen- 
sions agree well with exact solutions for adiabatic sound propagation and Couette 
flow with heat transfer. 
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The lattice Boltzmann (LB) method is a discrete, in space and time, microscopic, 
kinetic equation description for the evolution of the velocity distribution function of 
a fluid |T|, 0, ]!|. Like lattice gas (LG) automata [|J, LB methods are well suited for 
simulating a variety of physical systems in a parallel computing environment. As 
a result, the LB approach has found recent successes in a host of fluid dynamical 
problems, including flows in porous media [Q , magnetohydrodynamics || , immiscible 
fluids || and turbulence [J7], Its efficiency competes with, and in some cases exceeds, 
that of traditional numerical methods, while its physical interpretation is transparent. 

Noticeably absent, though, from the list of successful applications of LG and LB 
methods is a model which can simulate the full set of thermohydrodynamic equa- 
tions. Previous attempts at developing such a model have exclusively involved LG 
automata || whose Fermi-Dirac equilibrium distributions do not have sufficient 
flexibility to guarantee the correct form of the energy equation (3). LB methods are 
considerably more flexible, but have not, until now, been applied to this problem. 

The thermohydrodynamic equations of classical kinetic theory result from a Chapman- 
Enskog expansion of the continuum Boltzmann equation with the assumption of a 
Maxwellian equilibrium distribution. Since an exact Maxwellian distribution with a 
continuous distribution of velocities, both in angle and magnitude, cannot be imple- 
mented on a system that is discrete in both space and time, we seek an alternative 
distribution which will nevertheless give rise to the same macroscopic physics. In 
this Letter we address this issue and introduce a LB scheme which can simulate the 
following continuity, momentum, and energy equations for viscous, compressible, and 
heat-conducting flows: 
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where n is the fluid mass density, e is the internal energy per unit mass and is 
proportional to the temperature T, u is the local velocity, p is the pressure, and A, ji, 
and k are the second viscosity, shear viscosity, and thermal conductivity, respectively. 

The starting point of the LB method is the kinetic equation for the velocity dis- 
tribution function, /^(x, t): 

/ CTi (x + e ai ,t + 1) - / CTi (x,t) = Q ai , (4) 

where the nonnegative, real number / CT i(x, t) is the mass of fluid at lattice node x 
and time t, moving in direction i with speed, \e a i\ — cr — 1, 2, ...N, where N is the 
number of speeds. The a = speed corresponds to the component of the fluid which 
is at rest. The term Q ai represents the rate of change of f a i due to collisions. For 
computational efficiency, it is desirable to find the minimal set of a and i, for which 
a coarse-graining of the kinetic equation (4) leads to the macroscopic dynamics of 
interest . 

The microscopic dynamics associated with Equation (4) can be viewed as a two 
step process: free streaming and collision. During the free streaming step, f ai (yt-\-e ai ) 
is replaced by / CT j(x). Thus, each site exchanges mass with its neighbors, i.e. sites 
connected by lattice vectors e ai . In the collision step the distribution functions at 
each site then relax toward a state of local equilibrium. For simplicity, we consider the 



linearized, single-time-relaxation model of Bhatnagar, Gross, and Krook [|Tl|] , which 



has recently been applied to LB models p, 12, 13, 14 
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The collision operator Q ai conserves the local mass, momentum and kinetic energy: 
J2ai^ai = 0,Y.ai^aiG a i = 0, and ^1^1^/2 = 0, and the parameter, r, controls 
the rate at which the system relaxes to the local equilibrium, f^f. 

The LB method, unlike LGs, has considerable flexibility in the choice of the local 
equilibrium distribution. A general equilibrium distribution is given by a truncated 
power series in the local velocity u, valid for |u| <C 1, 

f% = A a + B a e ai ■ u + C a (e ai ■ u) 2 + D a u 2 + E a (e al • u) 3 + F a (e m ■ u)u 2 , (6) 

where the velocity is defined by: nu = ^aifai^ai- The coefficients, A, B, F, are 
functions of the local density n = Y,ai fai an d internal energy ne = J2ai fai( e ai — u ) 2 /2, 
and their functional forms depend on the geometry of the underlying lattice. 

The long-wavelength, low-frequency behavior of the this system is obtained by a 
Taylor series expansion of Equation (4) to second order in the lattice spacing and 
time step: 

dfai 1 d 1 d d 

— ^ + e ai ■ V/rt + -e ai e ai : Wf ai + e ai ■ V— f ai + -^r^rfvi = V ai . (7) 
at 2 at 2 at at 

In order to derive the macroscopic hydrodynamic equations, we adopt the following 
Chapman-Enskog multi-scale expansions. We expand the time derivative as 

9 d 9 d 

where the lower order terms in e vary more rapidly. Because we are interested in 
small departures from local equilibrium, we expand the distribution function as 

/« = /2 + e/j?+e 2 /£ ) + ..., (9) 



and the collision operator as 
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Substituting the above expansions into the kinetic equation, we find 



4rfS + eoi-Vf% = —f$ (11) 
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to order e, and 
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to order e 2 . With Equation (11) and some algebra, we can rewrite Equation (12) as: 

jfef% + (1 - ^)(^/- + e - • V ^) = ~flf. (13) 

Summing moments of Equations (11) and (13), we obtain to order e 2 , the continuity 
equation, 

^ + V-nu = 0, (14) 



the momemtum equation, 

and the energy equation, 

dne 



dt 



^ + v.n^o, (18) 



+ V • (neu) + V • q + P : Vu = 0. (16) 



The momentum flux tensor II = Y,ai[fai + (1 — /i^] e cri e <rij the heat flux, q Q 



(1/2) EcrAfai + (1 - ^)/il ) ]( e ai - u) 2 (e ai - u) a , and P is the pressure tensor, P Q/3 



(1/2) Eja + (1 " ^)/iP](e CT , - u) a (e CTi - uV 

To recover the Euler equations, we neglect the order e 2 terms and impose four fur- 
ther constraints on the equilibrium distribution function. The first of these constraints 
requires that the momentum flux tensor, 11^, be isotropic. The velocity independent 
portion of the tensor is then identified as the pressure, and this immediately results in 
the equation of state for an ideal gas, p = ne. The remaining two constraints require 
that the convective terms be Galilean invariant, and that the heat flux vanish to first 
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order in e, = 0. Thus we obtain the equations for compressible, inviscid and 
nonconducting flow of a monatomic gas. 

Retaining terms to order e 2 and imposing two additional constraints, we recover 
the Navier-Stokes level equations. These constraints are that the momentum flux 
tensor, n^ 1 ) , be isotropic and that the heat flux, , be proportional to the gradient 
of the temperature: VT. Note that the order e 2 terms describe diffusive 

processes, and, as assumed in Equation (8), evolve on a slower time scale than the 
convective terms associated with the order e Euler equations. 

To demonstrate the utility of the above LB method, we apply it to a two-dimensional 
triangular lattice. The model has one rest particle state, a — 0, for which e ai = 0, and 
two nonzero speeds for which e ai = o"(cos (2m/ 6), sin (2m/6)) for i — 0, 1 . . . , 6 and 
a = 1,2. The extension to three dimensions is straightforward and will be discussed 
elsewhere ||15|| . 

For this lattice geometry and the constraints discussed above, we can solve for the 
coefficients of the distribution function. One possible solution is the following: 

5 „ 2 a 4 4 2 . 1 2 1 

° = ~2 ne + n + 2ne ' Al = Q He ~ Q He ' A<2 = Q m ~ ^ n6 ' 

4 4 n 1 1 

B\ — -n ne, B 2 = -ne n, 

9 9 9 36 

m 8 4 ^ 11 
Ci = —n ne, Co = n H ne, 

9 3 ' 72 12 ' 

5 2 2 1 1 

D = — n + 2nt, Di = — n H — ne, D 2 = — n ne, 

4 ' 9 9 72 18 ' 

E 1 = —n, E 2 = —n, F t = 0, F 2 = 
27 108 

Identifying the coefficients in Equations (1) - (3) with the corresponding terms 
from the Chapman-Enskog expansion, we determine the values for the transport 
coefficients. The shear viscosity and the thermal conductivity are given by, \i = 
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ne(r — |) and k = 2ne(r — |), respectively, and yield a Prandtl number, Pr = 1/2. 
As in the case of a monatomic gas, the bulk viscosity vanishes because A = —fi. 

We carried out four numerical tests to determine the accuracy of this method for 
simulating Equations (l)-(3). Each test focused on one aspect of thermohydrody- 
namic transport. 

We determined the viscosity fi by simulating an isothermal Poiseuille flow. Numer- 
ical results demonstate that the model accurately reproduces a parabolic momentum 
profile (not shown here). The viscosity is related to the momentum at the channel 
center by fi = W 2 f/Su cen where / is the magnitude of the forcing, u cen is the velocity 



at the center, and W is the channel width [16|. In Figure 1, we show the dependence 



of viscosity on the relaxation time r and two internal energies e. The resulting vis- 
cosities from measurements agree with the Chapman- Enskog theory to within around 
1% over the entire range of parameters simulated. 

We determined the thermal conductivity k by measuring the heat transfer across a 
temperature gradient, using Fourier's law q = -kVT. By fixing the temperatures at 
the channel walls, we obtain a linear temperature profile and thus a constant gradient. 
Again, numerical results agree with the theoretical predictions quite well. Since the 
thermal conductivity has the same functional form as the viscosity, we also display 
the results in Figure 1. 

For the two dimensional LB scheme, linearized perturbation theory gives a simple 
relation between the adiabatic sound speed, c s and internal energy: c s = \/2e. In Fig- 
ure 2, we present the sound speed as a function of internal energy for both numerical 
measurements and theory - the agreement is evident. 

We simulated a Couette shear flow with a temperature gradient between the 
boundaries[|I7[]. For small temperature gradients the pressure is essentially constant 



across the channel, and the temperature profile has an analytic solution given by 



e* = (e — e )/(ei — e Q ) = (1/2) (1 + y*) + (Br/8)(1 — y* 2 ), where y* is the normalized 
distance from the center of the channel, t\ and eo are the internal energies of the 
upper and lower walls, respectively. The Brinkman number, Br is the product of 
the Prandtl and Eckert numbers. The agreement between theory and simulation, as 
shown in Figure 3, demonstrates the validity of the method in simulating flows in 
which energy dissipation is an important factor. 

In conclusion, we have developed a lattice Boltzmann scheme for the simulation of 
viscous, compressible, heat-conducting flows of an ideal monatomic gas. The kinetics 
of this model can be easily implemented on a parallel architecture machine. We 
have demonstrated theoretically and numerically that the macroscopic behavior of 
this model corresponds to that of Equations (1) - (3). Several issues remain. First, 
the current model uses the single-time-relaxation approximation, and this restricts 
simulations to flows with Prandtl number Pr = 1/2. In order to simulate flows with 
other Prandtl numbers, we should use a full matrix collision operator, which leads to 



a multi-time scale relaxation ||15|| . Second, the equation of state in the present model 
is that of ideal monatomic gas. To simulate non-ideal gases we may incorporate some 
internal degrees of freedom. An analysis of the numerical stability of the current 
model and its benchmarking against other computational fluid dynamical schemes 
are under investigation. 
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Figure Captions 

Fig.l : Shear viscosity (+) and thermal conductivity (O), as functions of relaxation 
parameter r. Upper curve corresponds to internal energy, e = 0.625, lower curve 
e = 0.5. The solid lines are the theoretical predictions. 

Fig. 2 : Numerical simulations of adiabatic sound speed (+) as a function of internal 
energy e. The solid line is the function \/2e. 

Fig. 3 : Normalized internal energy, e* for Couette flow with heat transfer for Brinkman 
numbers Br = 5 (+) and Br = 10 (□). The upper wall is moving with speed U\ = 0.1, 
and the lower wall is stationary. The solid lines represent analytical results. 
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